TP 2 - Extensions du modèle linéaire

Sélection de variables

Les données

On reprend le jeu de données “ozone.txt” du TP1.

Les 13 variables observées sont :

  • maxO3 : Maximum de concentration d’ozone observé sur la journée en \(\mu\)gr/m3
  • T9, T12, T15 : Température observée à 9, 12 et 15h
  • Ne9, Ne12, Ne15 : Nébulosité observée à 9, 12 et 15h
  • Vx9, Vx12, Vx15 : Composante E-O du vent à 9, 12 et 15h
  • maxO3v : Teneur maximum en ozone observée la veille
  • vent : orientation du vent à 12h
  • pluie : occurrence ou non de précipitations

Les données sont disponibles ici.

Pour les lire, vous pouvez utiliser la commande suivante :

Vous pouvez obtenir un résumé du jeu de données ozone avec la commande summary().

     maxO3              T9             T12             T15       
 Min.   : 42.00   Min.   :11.30   Min.   :14.00   Min.   :14.90  
 1st Qu.: 70.75   1st Qu.:16.20   1st Qu.:18.60   1st Qu.:19.27  
 Median : 81.50   Median :17.80   Median :20.55   Median :22.05  
 Mean   : 90.30   Mean   :18.36   Mean   :21.53   Mean   :22.63  
 3rd Qu.:106.00   3rd Qu.:19.93   3rd Qu.:23.55   3rd Qu.:25.40  
      Ne9             Ne12            Ne15           Vx9         
 Min.   :0.000   Min.   :0.000   Min.   :0.00   Min.   :-7.8785  
 1st Qu.:3.000   1st Qu.:4.000   1st Qu.:3.00   1st Qu.:-3.2765  
 Median :6.000   Median :5.000   Median :5.00   Median :-0.8660  
 Mean   :4.929   Mean   :5.018   Mean   :4.83   Mean   :-1.2143  
 3rd Qu.:7.000   3rd Qu.:7.000   3rd Qu.:7.00   3rd Qu.: 0.6946  
      Vx12             Vx15            maxO3v          vent      pluie   
 Min.   :-7.878   Min.   :-9.000   Min.   : 42.00   Est  :10   Pluie:43  
 1st Qu.:-3.565   1st Qu.:-3.939   1st Qu.: 71.00   Nord :31   Sec  :69  
 Median :-1.879   Median :-1.550   Median : 82.50   Ouest:50             
 Mean   :-1.611   Mean   :-1.691   Mean   : 90.57   Sud  :21             
 3rd Qu.: 0.000   3rd Qu.: 0.000   3rd Qu.:106.00                        
 [ getOption("max.print") est atteint -- ligne 1 omise ]

On réajuste ici un modèle de régression linéaire avec toutes les variables explicatives (cf TP1) pour la suite.


Call:
lm(formula = maxO3 ~ ., data = oz_quanti)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.566  -8.727  -0.403   7.599  39.458 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.24442   13.47190   0.909   0.3656    
T9          -0.01901    1.12515  -0.017   0.9866    
T12          2.22115    1.43294   1.550   0.1243    
T15          0.55853    1.14464   0.488   0.6266    
Ne9         -2.18909    0.93824  -2.333   0.0216 *  
Ne12        -0.42102    1.36766  -0.308   0.7588    
Ne15         0.18373    1.00279   0.183   0.8550    
Vx9          0.94791    0.91228   1.039   0.3013    
Vx12         0.03120    1.05523   0.030   0.9765    
Vx15         0.41859    0.91568   0.457   0.6486    
maxO3v       0.35198    0.06289   5.597 1.88e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.36 on 101 degrees of freedom
Multiple R-squared:  0.7638,    Adjusted R-squared:  0.7405 
F-statistic: 32.67 on 10 and 101 DF,  p-value: < 2.2e-16

Question 1

Enoncé

Dans le summary(reg.mul), un test est fait sur chaque coefficient. Il revient à tester que la variable n’apporte pas d’information supplémentaire sachant que toutes les autres variables sont dans le modèle.

Ici, nous souhaitons aller plus loin dans la possible simplification du modèle. Nous allons utiliser des procédures de choix de modèles pour sélectionner les variables significatives. On va ici comparer la sélection de variable obtenue par différents critères: BIC, AIC, \(R^2\) ajusté, \(C_p\) de Mallows. Pour cela, vous pouvez utiliser la fonction regsubsets() de la librairie leaps et la fonction stepAIC() de la librairie MASS. Commentez les résultats obtenus avec les différents critères, vous pourrez vous aider des commandes suivantes :

Correction

Start:  AIC=607.26
maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + 
    maxO3v

         Df Sum of Sq   RSS    AIC
- T9      1       0.1 20827 605.26
- Vx12    1       0.2 20827 605.26
- Ne15    1       6.9 20834 605.30
- Ne12    1      19.5 20847 605.36
- Vx15    1      43.1 20870 605.49
- T15     1      49.1 20876 605.52
- Vx9     1     222.6 21050 606.45
<none>                20827 607.26
- T12     1     495.5 21323 607.89
- Ne9     1    1122.6 21950 611.14
- maxO3v  1    6459.6 27287 635.51

Step:  AIC=605.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Vx12    1       0.1 20827 603.26
- Ne15    1       6.9 20834 603.30
- Ne12    1      22.1 20849 603.38
- Vx15    1      43.5 20871 603.49
- T15     1      49.4 20877 603.52
- Vx9     1     264.0 21091 604.67
<none>                20827 605.26
- T12     1     641.4 21469 606.66
- Ne9     1    1183.6 22011 609.45
- maxO3v  1    7112.2 27940 636.16

Step:  AIC=603.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Ne15    1       6.9 20834 601.30
- Ne12    1      23.1 20850 601.38
- T15     1      50.0 20877 601.53
- Vx15    1      79.6 20907 601.69
- Vx9     1     320.1 21148 602.97
<none>                20827 603.26
- T12     1     647.3 21475 604.69
- Ne9     1    1185.1 22013 607.46
- maxO3v  1    7112.6 27940 634.16

Step:  AIC=601.3
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Ne12    1      16.2 20851 599.38
- T15     1      45.2 20880 599.54
- Vx15    1      75.1 20910 599.70
- Vx9     1     325.2 21160 601.03
<none>                20834 601.30
- T12     1     971.8 21806 604.40
- Ne9     1    1215.8 22050 605.65
- maxO3v  1    7106.3 27941 632.17

Step:  AIC=599.38
maxO3 ~ T12 + T15 + Ne9 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- T15     1      45.6 20896 597.63
- Vx15    1      77.6 20928 597.80
- Vx9     1     337.6 21188 599.18
<none>                20851 599.38
- T12     1    1034.0 21885 602.80
- Ne9     1    2155.5 23006 608.40
- maxO3v  1    7130.2 27981 630.33

Step:  AIC=597.63
maxO3 ~ T12 + Ne9 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Vx15    1      74.1 20970 596.02
- Vx9     1     366.5 21263 597.58
<none>                20896 597.63
- Ne9     1    2241.1 23137 607.04
- T12     1    6712.6 27609 626.83
- maxO3v  1    7381.5 28278 629.51

Step:  AIC=596.02
maxO3 ~ T12 + Ne9 + Vx9 + maxO3v

         Df Sum of Sq   RSS    AIC
<none>                20970 596.02
- Vx9     1     903.4 21874 598.75
- Ne9     1    2714.8 23685 607.66
- T12     1    6650.4 27621 624.88
- maxO3v  1    7363.5 28334 627.73

Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = oz_quanti)

Coefficients:
(Intercept)          T12          Ne9          Vx9       maxO3v  
    12.6313       2.7641      -2.5154       1.2929       0.3548  

Avec le \(C_p\) de Mallows, le critère BIC et le R2 ajusté, on retient 4 variables explicatives : T12, Ne9, Vx9 et maxO3v.
Même concluson en utilisant la fonction stepAIC.

On retrouve une cohérence avec les corrélations entre variables

Question 2

Enoncé

A la question précédente, nous retenons 4 variables explicatives T12, Ne9, Vx9 et maxO3v.
A l’aide des commandes suivantes, testez le sous-modèle avec les 4 variables retenues contre le modèle complet. Commentez.

Correction


Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = oz_quanti)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.396  -8.377  -1.086   7.951  40.933 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.63131   11.00088   1.148 0.253443    
T12          2.76409    0.47450   5.825 6.07e-08 ***
Ne9         -2.51540    0.67585  -3.722 0.000317 ***
Vx9          1.29286    0.60218   2.147 0.034055 *  
maxO3v       0.35483    0.05789   6.130 1.50e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14 on 107 degrees of freedom
Multiple R-squared:  0.7622,    Adjusted R-squared:  0.7533 
F-statistic: 85.75 on 4 and 107 DF,  p-value: < 2.2e-16
Analysis of Variance Table

Model 1: maxO3 ~ T12 + Ne9 + Vx9 + maxO3v
Model 2: maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + 
    maxO3v
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    107 20970                           
2    101 20827  6    143.01 0.1156 0.9944

La p-valeur est \(p=0.994 \geq 0.05\) donc on ne rejette pas l’hypothèse nulle au risque \(5\%\) donc on retient le modèle simplifié avec 4 variables explicatives.

Régressions régularisées

On va s’intéresser aux régressions ridge, Lasso et Elastic-net dans cette partie. On commence par centrer et réduire les données avant de mettre en place et comparer ces méthodes de régression régularisées.

Régression Ridge

Question 1

Enoncé

A l’aide de la fonction glmnet(), ajustez une régression ridge en faisant varier \(\tau\) sur une grille. On stockera le résultat dans la variable fitridge. Explorez le contenu de fitridge.

Correction

          Length Class     Mode   
a0         1001  -none-    numeric
beta      10010  dgCMatrix S4     
df         1001  -none-    numeric
dim           2  -none-    numeric
lambda     1001  -none-    numeric
dev.ratio  1001  -none-    numeric
nulldev       1  -none-    numeric
npasses       1  -none-    numeric
jerr          1  -none-    numeric
offset        1  -none-    logical
call          7  -none-    call   
nobs          1  -none-    numeric

fitridge$lambda contient les valeurs du paramètre \(\tau\) dans l’ordre décroissant. fitridge$beta contient la valeur des paramètres \(\theta_j\) pour chaque valeur de \(\tau\). Comme on a centré les données, on n’a pas d’intercept donc fitridge$a0 est un vecteur de zéro.

Pour l’interprétation des autres sorties, vous pouvez vous reporter à l’aide de la fonction glmnet.

Question 2

Correction

Question 3

Régression Lasso

Question 1

Enoncé

A l’aide de la fonction glmnet(), ajustez une régression Lasso en faisant varier \(\tau\) sur une grille. On stockera le résultat dans la variable fitlasso. Explorez le contenu de fitlasso.

Correction

          Length Class     Mode   
a0         1001  -none-    numeric
beta      10010  dgCMatrix S4     
df         1001  -none-    numeric
dim           2  -none-    numeric
lambda     1001  -none-    numeric
dev.ratio  1001  -none-    numeric
nulldev       1  -none-    numeric
npasses       1  -none-    numeric
jerr          1  -none-    numeric
offset        1  -none-    logical
call          7  -none-    call   
nobs          1  -none-    numeric

Voici quelques exemples de valeurs estimées pour les paramètres pour plusieurs valeurs de \(\tau\)

Question 2

Correction

Question 3

Enoncé

A l’aide de la fonction cv.glmnet() mettez en place une validation croisée pour sélectionner le “meilleur” \(\tau\) par MSE.

La valeur de \(\tau\) sélectionnée vaut ….

Question 4

Enoncé

Quelle sélection de variables obtient-on alors ?

Correction

11 x 1 sparse Matrix of class "dgCMatrix"
                      s1
(Intercept)  .          
T9           .          
T12          0.300068858
T15          0.079165203
Ne9         -0.144351908
Ne12        -0.002279795
Ne15         .          
Vx9          0.039391227
Vx12         .          
Vx15         .          
maxO3v       0.251543269

Régression Elastic Net

Question 1

Correction

Question 2

Enoncé

Comparez les coefficients obtenus avec la régression linéaire, la régression ridge, la régression Lasso et la régression ElasticNet

Correction

Modèle linéaire généralisé

Régression logistique

Dans cette section, nous allons nous intéresser au jeu de données SAheart disponible dans la librairie bestglm. On souhaite expliquer la présence/absence d’une maladie cardiovasculaire (chd) en fonction de 9 variables et on dispose pour cela d’un échantillon de \(n=462\) individus.

Question 1

A l’aide des commandes suivantes, chargez les données et faites quelques statistiques descriptives pour vous familiariser avec les données.

La description des données est disponible dans l’aide de R (tapez ?SAheart dans la console).

Le jeu de données est composé des variables suivantes :

  1. sbp: pression artérielle systolique
  2. tabac: tabac cumulatif (kg)
  3. ldl: cholestérol des lipoprotéines de basse densité
  4. adiposité:
  5. famhist: antécédents familiaux de maladie cardiaque (Présent, Absent)
  6. type: comportement de type-A
  7. obésité:
  8. alcool: consommation actuelle d’alcool
  9. âge: âge au début
  10. chd: réponse, maladie coronarienne
'data.frame':   462 obs. of  10 variables:
 $ sbp      : int  160 144 118 170 134 132 142 114 114 132 ...
 $ tobacco  : num  12 0.01 0.08 7.5 13.6 6.2 4.05 4.08 0 0 ...
 $ ldl      : num  5.73 4.41 3.48 6.41 3.5 6.47 3.38 4.59 3.83 5.8 ...
 $ adiposity: num  23.1 28.6 32.3 38 27.8 ...
 $ famhist  : Factor w/ 2 levels "Absent","Present": 2 1 2 2 2 2 1 2 2 2 ...
 $ typea    : int  49 55 52 51 60 62 59 62 49 69 ...
 $ obesity  : num  25.3 28.9 29.1 32 26 ...
 $ alcohol  : num  97.2 2.06 3.81 24.26 57.34 ...
 $ age      : int  52 63 46 58 49 45 38 58 29 53 ...
 $ chd      : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 1 2 1 2 ...
      sbp           tobacco             ldl           adiposity    
 Min.   :101.0   Min.   : 0.0000   Min.   : 0.980   Min.   : 6.74  
 1st Qu.:124.0   1st Qu.: 0.0525   1st Qu.: 3.283   1st Qu.:19.77  
 Median :134.0   Median : 2.0000   Median : 4.340   Median :26.11  
 Mean   :138.3   Mean   : 3.6356   Mean   : 4.740   Mean   :25.41  
 3rd Qu.:148.0   3rd Qu.: 5.5000   3rd Qu.: 5.790   3rd Qu.:31.23  
 Max.   :218.0   Max.   :31.2000   Max.   :15.330   Max.   :42.49  
    famhist        typea         obesity         alcohol            age       
 Absent :270   Min.   :13.0   Min.   :14.70   Min.   :  0.00   Min.   :15.00  
 Present:192   1st Qu.:47.0   1st Qu.:22.98   1st Qu.:  0.51   1st Qu.:31.00  
               Median :53.0   Median :25.80   Median :  7.51   Median :45.00  
               Mean   :53.1   Mean   :26.04   Mean   : 17.04   Mean   :42.82  
               3rd Qu.:60.0   3rd Qu.:28.50   3rd Qu.: 23.89   3rd Qu.:55.00  
               Max.   :78.0   Max.   :46.58   Max.   :147.19   Max.   :64.00  
 chd    
 0:302  
 1:160  
        
        
        
        

Question 2

Enoncé

Ajustez un modèle de régression logistique multiple additif à l’aide de la fonction glm().

Correction

Modèle de régression logistique :

\[chd_i\sim\mathcal{B}(\pi_i)\] avec
\[g(\pi_i)=\ln\left(\frac{\pi_i}{1-\pi_i}\right) = \theta_0 + \sum_{j=1}^{9} \theta_j x_i^{(j)} \textrm{ dont } x_i^{(famhist)} = \mathbb{1}_{famhist_i = present}\]


Call:
glm(formula = chd ~ ., family = binomial(link = "logit"), data = SAheart)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7781  -0.8213  -0.4387   0.8889   2.5435  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -6.1507209  1.3082600  -4.701 2.58e-06 ***
sbp             0.0065040  0.0057304   1.135 0.256374    
tobacco         0.0793764  0.0266028   2.984 0.002847 ** 
ldl             0.1739239  0.0596617   2.915 0.003555 ** 
adiposity       0.0185866  0.0292894   0.635 0.525700    
famhistPresent  0.9253704  0.2278940   4.061 4.90e-05 ***
typea           0.0395950  0.0123202   3.214 0.001310 ** 
obesity        -0.0629099  0.0442477  -1.422 0.155095    
alcohol         0.0001217  0.0044832   0.027 0.978350    
age             0.0452253  0.0121298   3.728 0.000193 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 596.11  on 461  degrees of freedom
Residual deviance: 472.14  on 452  degrees of freedom
AIC: 492.14

Number of Fisher Scoring iterations: 5

Question 3

Enoncé

Déterminez la valeur du pseudo-\(R^2\) pour le modèle modellogit.

Correction

[1] 0.2079628

La valeur du pseudo-\(R^2\) n’est pas très élevé.

Question 4

Enoncé

Commentez les résultats obtenus via la commande suivante :

Correction

                    Estimate  Std. Error     z value     Pr(>|z|)
(Intercept)    -6.1507208650 1.308260018 -4.70145138 2.583188e-06
sbp             0.0065040171 0.005730398  1.13500273 2.563742e-01
tobacco         0.0793764457 0.026602843  2.98375801 2.847319e-03
ldl             0.1739238981 0.059661738  2.91516648 3.554989e-03
adiposity       0.0185865682 0.029289409  0.63458325 5.257003e-01
famhistPresent  0.9253704194 0.227894010  4.06052980 4.896149e-05
typea           0.0395950250 0.012320227  3.21382267 1.309805e-03
obesity        -0.0629098693 0.044247743 -1.42176449 1.550946e-01
alcohol         0.0001216624 0.004483218  0.02713729 9.783502e-01
age             0.0452253496 0.012129752  3.72846442 1.926501e-04

On a en dernière colonne les pvaleurs associées à chaque \(Z\)-test pour tester la nullité de chacun des coefficients. On peut remarquer que les variables sbp, adiposity, obesity et alcohol ont une pvaleur assez élevée. On va poursuivre à l’étude en cherchant à simplifier le modèle par sélection de variables.

Question 5

Correction

On commence par une procédure de sélection de variables avec la fonction bestglm() en utilisant le critère BIC :


Call:  glm(formula = y ~ ., family = family, data = Xi, weights = weights)

Coefficients:
   (Intercept)         tobacco             ldl  famhistPresent           typea  
      -6.44644         0.08038         0.16199         0.90818         0.03712  
           age  
       0.05046  

Degrees of Freedom: 461 Total (i.e. Null);  456 Residual
Null Deviance:      596.1 
Residual Deviance: 475.7    AIC: 487.7
Fitting algorithm:  BIC-glm
Best Model:
            df deviance
Null Model 456 475.6856
Full Model 461 596.1084

    likelihood-ratio test - GLM

data:  H0: Null Model vs. H1: Best Fit BIC-glm
X = 120.42, df = 5, p-value < 2.2e-16

Le modèle retenu est \[ logit(\pi_i) = \theta_0 + \theta_2\ tobacco_i + \theta_3\ ldl_i + \theta_5\ \mathbb{1}_{famhist_i = present} + \theta_6\ typea_i + \theta_9\ age_i \]

On fait de même avec le critère AIC :


Call:  glm(formula = y ~ ., family = family, data = Xi, weights = weights)

Coefficients:
   (Intercept)         tobacco             ldl  famhistPresent           typea  
      -6.44644         0.08038         0.16199         0.90818         0.03712  
           age  
       0.05046  

Degrees of Freedom: 461 Total (i.e. Null);  456 Residual
Null Deviance:      596.1 
Residual Deviance: 475.7    AIC: 487.7

On sélectionne le même modèle.

On peut aussi faire la procédure de sélection de variables à l’aide de la fonction step(). Par défaut, on utilise le critère AIC. On peut voir les étapes de la procédure de sélection backward. On retrouve le même sous-modèle retenu.

Start:  AIC=492.14
chd ~ sbp + tobacco + ldl + adiposity + famhist + typea + obesity + 
    alcohol + age

            Df Deviance    AIC
- alcohol    1   472.14 490.14
- adiposity  1   472.55 490.55
- sbp        1   473.44 491.44
<none>           472.14 492.14
- obesity    1   474.23 492.23
- ldl        1   481.07 499.07
- tobacco    1   481.67 499.67
- typea      1   483.05 501.05
- age        1   486.53 504.53
- famhist    1   488.89 506.89

Step:  AIC=490.14
chd ~ sbp + tobacco + ldl + adiposity + famhist + typea + obesity + 
    age

            Df Deviance    AIC
- adiposity  1   472.55 488.55
- sbp        1   473.47 489.47
<none>           472.14 490.14
- obesity    1   474.24 490.24
- ldl        1   481.15 497.15
- tobacco    1   482.06 498.06
- typea      1   483.06 499.06
- age        1   486.64 502.64
- famhist    1   488.99 504.99

Step:  AIC=488.55
chd ~ sbp + tobacco + ldl + famhist + typea + obesity + age

          Df Deviance    AIC
- sbp      1   473.98 487.98
<none>         472.55 488.55
- obesity  1   474.65 488.65
- tobacco  1   482.54 496.54
- ldl      1   482.95 496.95
- typea    1   483.19 497.19
- famhist  1   489.38 503.38
- age      1   495.48 509.48

Step:  AIC=487.98
chd ~ tobacco + ldl + famhist + typea + obesity + age

          Df Deviance    AIC
- obesity  1   475.69 487.69
<none>         473.98 487.98
- tobacco  1   484.18 496.18
- typea    1   484.30 496.30
- ldl      1   484.53 496.53
- famhist  1   490.58 502.58
- age      1   502.11 514.11

Step:  AIC=487.69
chd ~ tobacco + ldl + famhist + typea + age

          Df Deviance    AIC
<none>         475.69 487.69
- ldl      1   484.71 494.71
- typea    1   485.44 495.44
- tobacco  1   486.03 496.03
- famhist  1   492.09 502.09
- age      1   502.38 512.38

Pour utiliser le critère BIC avec la fonction step(), il faut utiliser la syntaxe suivante :

Question 6

Enoncé

Confirmez par un test de sous-modèle le sous-modèle retenu modelbest dans la question précédente.

Par quels critères peut-on aussi comparer les deux modèles modelbest et modellogit ?

Correction

On commence par ajuster le modèle retenu par les méthodes de sélection de variables


Call:
glm(formula = chd ~ tobacco + ldl + famhist + typea + age, family = binomial, 
    data = SAheart)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9165  -0.8054  -0.4430   0.9329   2.6139  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -6.44644    0.92087  -7.000 2.55e-12 ***
tobacco         0.08038    0.02588   3.106  0.00190 ** 
ldl             0.16199    0.05497   2.947  0.00321 ** 
famhistPresent  0.90818    0.22576   4.023 5.75e-05 ***
typea           0.03712    0.01217   3.051  0.00228 ** 
age             0.05046    0.01021   4.944 7.65e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 596.11  on 461  degrees of freedom
Residual deviance: 475.69  on 456  degrees of freedom
AIC: 487.69

Number of Fisher Scoring iterations: 5

On peut par exemple remarquer que le test de nullité de chacun des coefficients séparement a une p-valeur \(\ll 0.05\).

On fait ensuite un test de sous-modèle entre modelbest et modellogit

Analysis of Deviance Table

Model 1: chd ~ tobacco + ldl + famhist + typea + age
Model 2: chd ~ sbp + tobacco + ldl + adiposity + famhist + typea + obesity + 
    alcohol + age
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       456     475.69                     
2       452     472.14  4   3.5455    0.471

La pvaleur valant 0.47, on ne rejette pas le sous-modèle au risque \(5\%\). Pour le sous-modèle on a \(k_0=6\) paramètres à estimer donc \(n-k_0=462-6=456\)
Pour le modèle complet on a \(k_1=10\) paramètres à estimer donc \(n-k_0=462-10=452\)
\(\mathcal{D}(M_0) - \mathcal{D}(M_1) \stackrel{\mathcal{L}}{\longrightarrow} \chi^2(k_1-k_0 =4)\)

Critère pseudo-\(R^2\) :

[1] 0.202015
[1] 0.2079628

Critère AIC :

[1] 487.6856
[1] 492.14

Question 7

Enoncé

Comment s’interprètent les différents coefficients de modelbest ?

   (Intercept)        tobacco            ldl famhistPresent          typea 
   0.001586152    1.083693731    1.175850406    2.479793436    1.037812583 
           age 
   1.051755195 

Correction

Pour l’interprétation des paramètres, il faut considérer des odds-ratio, attention entre variables qualitatives et quantitatives.

Régression loglinéaire

Dans cette partie, on souhaite étudier la diversité des fourmis sur le site expérimental des Nourragues en Guyane Française. Les données, disponibles dans le fichier Fourmis.txt, sont issues du protocole expérimental suivant : des morceaux de 1\(m^2\) de litière ont été récoltés. Ils ont été pesés (car le poids de litière est vu comme un indicateur de l’épaisseur de la litière) et le nombre d’espèces différentes présentes dans l’échantillon a été dénombré. Enfin les conditions de recueil (humides ou sèches) ont été notées pour tester leur influence sur la présence des fourmis. L’objectif est de mettre en évidence les variables qui influencent potentiellement le nombre d’espèces de fourmis présentes dans le milieu.

Question 1

Chargez le jeu de données Fourmis et executez les commandes suivantes pour faire quelques statistiques descritpives.

   Effectifs         Weight         Site    Conditions
 Min.   : 2.00   Min.   : 6.200   FLWT:50   Dry:83    
 1st Qu.: 7.00   1st Qu.: 8.525   FTWT:50   Wet:87    
 Median :10.00   Median : 9.950   GPWT:50             
 Mean   :11.19   Mean   :10.168   INWT:20             
 3rd Qu.:15.00   3rd Qu.:11.600                       
 Max.   :28.00   Max.   :17.800                       
# A tibble: 8 x 3
# Groups:   Site [4]
  Site  Conditions `n()`
  <fct> <fct>      <int>
1 FLWT  Dry           24
2 FLWT  Wet           26
3 FTWT  Dry           24
4 FTWT  Wet           26
5 GPWT  Dry           22
6 GPWT  Wet           28
7 INWT  Dry           13
8 INWT  Wet            7

Question 2

Enoncé

Dans cette question, on cherche à expliquer le nombre de fourmis présentes dans un échantillon de sol en prenant en compte les conditions de recueil, le site et l’interaction entre ces deux facteurs.

  • Ecrivez un modèle linéaire généralisé modpois adapté et programmez-le à l’aide de la fonction glm().
  • Etudiez si on peut simplifier le modèle modpois. Vous pouvez vous aider de la fonction anova(....,test="Chisq").

  • A l’aide des commandes suivantes, on obtient le nombre moyen d’espèces de fourmis attendu sur un échantillon de terre aux différents sites pour les deux conditions. Commentez.

Correction

On considère un modèle loglinéaire modpois pour expliquer le nombre de fourmis présentes dans un échantillon de sol en prenant en compte les conditions de recueil, le site et l’interaction entre ces deux facteurs.

\[ Eff_i\sim\mathcal{P}(\lambda_i) \textrm{ et } \ln(\lambda_i) = \theta_0 + \sum_{j\in S}(\theta_j + \gamma_j \mathbb{1}_{Conditions_i=WET})\mathbb{1}_{Site_i=j} + \alpha \mathbb{1}_{Conditions_i=WET} \] avec \(S=\{FTWT,GPWT,INTW\}\).


Call:
glm(formula = Effectifs ~ Site * Conditions, family = poisson, 
    data = Fourmis)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.7219  -1.2101  -0.2415   0.9066   2.9152  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)             2.82583    0.04969  56.869  < 2e-16 ***
SiteFTWT               -0.62861    0.08425  -7.461 8.60e-14 ***
SiteGPWT               -0.69113    0.08857  -7.803 6.06e-15 ***
SiteINWT               -0.42794    0.09727  -4.399 1.09e-05 ***
ConditionsWet          -0.13851    0.07132  -1.942   0.0521 .  
SiteFTWT:ConditionsWet -0.17515    0.12476  -1.404   0.1603    
SiteGPWT:ConditionsWet  0.28473    0.11880   2.397   0.0165 *  
SiteINWT:ConditionsWet  0.63099    0.14148   4.460 8.20e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 530.36  on 169  degrees of freedom
Residual deviance: 329.68  on 162  degrees of freedom
AIC: 1045.4

Number of Fisher Scoring iterations: 4

On étudie si on peut simplifier le modèle modpois. On commence par tester si on peut enlever l’interaction.

Analysis of Deviance Table

Model 1: Effectifs ~ Site + Conditions
Model 2: Effectifs ~ Site * Conditions
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       165     361.05                          
2       162     329.68  3   31.367 7.114e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On rejette l’hypothèse nulle, on conserve l’effet d’interaction. On ne peut donc pas simplifier le modèle.

Nombre moyen d’espèces de fourmis attendu sur un échantillon de terre aux différents sites pour les deux conditions :

  Site Conditions lambdahat
1 FLWT        Dry 16.875000
2 FTWT        Dry  9.000000
3 GPWT        Dry  8.454545
4 INWT        Dry 11.000000
5 FLWT        Wet 14.692308
6 FTWT        Wet  6.576923
7 GPWT        Wet  9.785714
8 INWT        Wet 18.000000

Question 3

Enoncé

On souhaite maintenant étudier le nombre d’espèces de fourmis présentes dans une unité de volume, en fonction des caractéristiques du site. Nous allons ici intégrer l’information du poids (Weight) qui est vu comme un indicateur de l’épaisseur de la litière.

On remarque que la relation entre les paramètres des lois de Poisson du nombre d’espèces présentes en un site et du nombre d’espèces présentes dans une unité de volume est

\[ \mathbb{E}[\textrm{nb esp.}]=\mathbb{E}[\textrm{nb esp. par unité de volume}]\times Weight \]

On en déduit donc un modèle pour le nombre d’espèces par unité de volume :

\[Eff_i\sim\mathcal{P}(\lambda_i)\] avec

\[\ln(\lambda_i) = \ln(Weight_i) + \theta_0 + \sum_{j\in S}(\theta_j + \gamma_j \mathbb{1}_{Conditions_i=WET})\mathbb{1}_{Site_i=j} + \alpha \mathbb{1}_{Conditions_i=WET} \]

La variable Weight est donc considérée comme un offset.

  • Ajustez ce modèle à l’aide de l’option "offset" de la fonction glm().

  • Quel est le nombre moyen d’espèces de fourmis attendu sur un échantillon de terre qui pèse 10kg aux différents sites pour les deux conditions ?

Correction

Modèle pour le nombre d’espèces par unité de volume : \[Eff_i\sim\mathcal{P}(\lambda_i)\] avec

\[\ln(\lambda_i) = \ln(Weight_i) + \theta_0 + \sum_{j\in S}(\theta_j + \gamma_j \mathbb{1}_{Conditions_i=WET})\mathbb{1}_{Site_i=j} + \alpha \mathbb{1}_{Conditions_i=WET} \]


Call:
glm(formula = Effectifs ~ Site * Conditions, family = "poisson", 
    data = Fourmis, offset = log(Weight))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2267  -0.9268  -0.1522   0.8252   3.6257  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)             0.49531    0.04969   9.968  < 2e-16 ***
SiteFTWT               -0.55200    0.08425  -6.552 5.69e-11 ***
SiteGPWT               -0.74225    0.08857  -8.380  < 2e-16 ***
SiteINWT               -0.38293    0.09727  -3.937 8.26e-05 ***
ConditionsWet          -0.02802    0.07132  -0.393  0.69438    
SiteFTWT:ConditionsWet -0.37426    0.12476  -3.000  0.00270 ** 
SiteGPWT:ConditionsWet  0.16746    0.11880   1.410  0.15866    
SiteINWT:ConditionsWet  0.47387    0.14148   3.349  0.00081 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 534.13  on 169  degrees of freedom
Residual deviance: 301.41  on 162  degrees of freedom
AIC: 1017.2

Number of Fisher Scoring iterations: 4

Le pseudo-R\(^2\) de ce modèle :

[1] 0.4357083

Le nombre moyen d’espèces de fourmis attendu sur un échantillon de terre qui pèse 10kg aux différents sites pour les deux conditions ?

  Site Conditions Weight lambdahat
1 FLWT        Dry     10 16.410049
2 FTWT        Dry     10  9.448819
3 GPWT        Dry     10  7.811844
4 INWT        Dry     10 11.189358
5 FLWT        Wet     10 15.956558
6 FTWT        Wet     10  6.319290
7 GPWT        Wet     10  8.980662
8 INWT        Wet     10 17.475728

Question 4

Enoncé

Est-ce-qu’une modélisation avec un modèle prenant en compte la sur-dispersion est plus appropriée ?

Correction


Call:
glm(formula = Effectifs ~ Site * Conditions, family = quasipoisson(link = "log"), 
    data = Fourmis, offset = log(Weight))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2267  -0.9268  -0.1522   0.8252   3.6257  

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.49531    0.06739   7.350 9.23e-12 ***
SiteFTWT               -0.55200    0.11426  -4.831 3.12e-06 ***
SiteGPWT               -0.74225    0.12012  -6.179 5.02e-09 ***
SiteINWT               -0.38293    0.13191  -2.903  0.00421 ** 
ConditionsWet          -0.02802    0.09672  -0.290  0.77239    
SiteFTWT:ConditionsWet -0.37426    0.16918  -2.212  0.02836 *  
SiteGPWT:ConditionsWet  0.16746    0.16110   1.039  0.30015    
SiteINWT:ConditionsWet  0.47387    0.19186   2.470  0.01455 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 1.839056)

    Null deviance: 534.13  on 169  degrees of freedom
Residual deviance: 301.41  on 162  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Pseudo-R\(^2\) pour ce modèle :

[1] 0.4357083